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Abstract. Harmonic inversion has already been proven to be a powerful tool for the analysis of quantum 
spectra and the periodic orbit orbit quantization of chaotic systems. The harmonic inversion technique 
circumvents the convergence problems of the periodic orbit sum and the uncertainty principle of the usual 
Fourier analysis, thus yielding results of high resolution and high precision. Based on the close analogy 
between periodic orbit trace formulae for regular and chaotic systems the technique is generalized in this 
paper for the semiclassical quantization of integrable systems. Thus, harmonic inversion is shown to be a 
universal tool which can be applied to a wide range of physical systems. The method is further generalized in 
two directions: Firstly, the periodic orbit quantization will be extended to include higher order h corrections 
to the periodic orbit sum. Secondly, the use of cross-correlated periodic orbit sums allows us to significantly 
reduce the required number of orbits for semiclassical quantization, i.e., to improve the efficiency of the 
semiclassical method. As a representative of regular systems, we choose the circle billiard, whose periodic 
orbits and quantum eigenvalues can easily be obtained. 

PACS. 03.65.Sq Semiclassical theories and applications 



1 Introduction 

A question of fundamental interest for systems with both 
regular and chaotic dynamics is how quantum mechanical 
eigenvalues can be obtained by quantization of classical 
orbits. The EBK torus quantization method of Einstein, 
Brillouin, and Keller fl],§,|| is restricted to integrable sys- 
tems, i.e., the method cannot be generalized to systems 
with an underlying chaotic or mixed regular-chaotic dy- 
namics Q. Furthermore, EBK quantization requires the 
knowledge of all the constants of motion, which are not 
normally given in explicit form, and therefore practical 
EBK quantization based on the direct or indirect numer- 
ical construction of the constants of motion turns out to 
be a formidable task As an alternative, EBK quan- 
tization was recast as a sum over all periodic orbits of a 
given topology on respective tori by Berry and Tabor Q . 
In contrast to EBK-quantization, periodic orbit theory can 
be applied to systems with more general classical dynam- 
ics: Gutzwiller's trace formula [^,|^| for chaotic systems 
and the corresponding Berry- Tabor formula for regular 
systems || give the semiclassical approximation for the 
density of states as a sum over the periodic orbits of the 
underlying classical system. However, a fundamental prob- 
lem of these periodic orbit sums is that they usually do 
not converge, or if they do, the convergence is extremely 
slow. During recent years, various techniques have been 
developed to overcome this problem. Most of them are 
especially designed for chaotic systems MM, [L0| and can- 



not be applied to systems with regular or mixed regular- 
chaotic dynamics, or they depend on special properties of 
the system such as the existence of a symbolic dynamics. 
They are therefore restricted to a relatively small num- 
ber of physical systems. It would be desirable to have a 
method at hand which is universal in the sense that it is 
applicable for all types of underlying classical dynamics. 

Recently, a method for periodic orbit quantization, 
based on harmonic inversion of a semiclassical signal 
has been developed and successfully applied to classically 
chaotic systems Jll], 12],[l]| . The aim of the present paper 
is to demonstrate that this technique is equally power- 
ful in reproducing the spectra of regular systems. The 
semiclassical quantization of integrable and chaotic sys- 
tems on an equal footing will be the basis for applications 
to systems with even more general, i.e., mixed regular- 
chaotic dynamics [Q. Furthermore, the harmonic inver- 
sion technique is generalized in two directions: Firstly, the 
periodic orbit quantization will be extended to include 
higher order h corrections [15]], and, second ly, the use of 
cross-correlated periodic or 
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allows us 

to significantly reduce the required number of orbits for 
semiclassical quantization, i.e., to improve the efficiency 
of the semiclassical method. As a representative of regu- 
lar systems, we choose the circle billiard whose periodic 
orbits and quantum eigenvalues can easily be obtained. 
The paper is organized as follows: 

In Section || we give a brief overview over the periodic 
orbit theory for integrable systems, especially the Berry- 



2 Kirsten Weibert et al.: Use of harmonic inversion techniques in the periodic orbit quantization of integrable systems 



Tabor formula, which is the analogue for integrable sys- 
tems to Gutzwiller's trace formula for chaotic systems. We 
then calculate the explicit expression for the density of 
states of the circle billiard from the Berry- Tabor formula. 
The equations are generalized in two directions, firstly, to 
the density of states weighted with the diagonal matrix el- 
ements of one or more given operators Jl6| , and, secondly, 
to include higher order h corrections in the periodic orbit 
sum |p^| . 

The high precision analysis of quantum spectra and the 
method for the analytic continuation of non-convergent 
periodic orbit sums applied in this paper are based 
on the harmonic inversion of time signals. In Sec- 
tion m we briefly introduce harmonic inversion by filter- 
diagonalization We also discuss an extension of 

the filter-diagonalization method to cross-correlation func- 
tions |n],|2l],|22|], which can be used to extract semiclassi- 
cal eigenvalues and matrix elements from cross-correlated 
periodic orbit sums with a significantly reduced set of pe- 
riodic orbits fl7| . 

Harmonic inversion circumvents the uncertainty prin- 
ciple of the conventional Fourier transform and can be 
used for the high precision analysis of quantum spectra 
||l3|,^3). In Section |J the method will be applied to the 
quantum spectra of the circle billiard. The analysis will 
verify the validity of the Berry- Tabor formula and its gen- 
eralization to spectra weight ed with diagonal matrix el- 
ements discussed in Section 2.E. Furthermore, harmonic 



inversion will be applied to determine the higher order h 
contributions to the periodic orbit sum. The Gutzwiller 
and the Berry- Tabor formula are only the leading or- 
der contributions of an expansion of the density of states 
in terms of fi and therefore only yield semiclassical ap- 
proximations to the eigenvalues. By analyzing the differ- 
ence spectrum between exact and semiclassical eigenval- 
ues, first order H corrections to the periodic orbit sum can 
be determined, as we will demonstrate in Section 4.2. The 



results are compared with the analytic expressions for the 
h expansion of the periodic orbit sum given in Section 2A . 

In Section ||, we turn to the periodic orbit quantiza- 
tion of integrable systems. Firstly, we show how in general 
the problem of extracting semiclassical eigenvalues from 
periodic orbit sums can be reformulated as a harmonic 
inversion problem: A semiclassical signal is constructed 
from the periodic orbit sum, the analysis of which yields 
the semiclassical eigenvalues of the system. The general 
formulae are then applied to the circle billiard and the 
results are compared t o th e exact quantum and the EBK 
eigenvalues. In Section 5.3 it is demonstrated how the ac- 
curacy of the semiclassical eigenvalues can be significantly 
improved with the help of higher order h corrections to the 
periodic orbit sum. 



In Section 5.4 we address the question of how to 
improve the efficiency of the semiclassical quantization 
method, i.e., how to extract the same number of eigenval- 
ues with a reduced set of periodic orbits, which is impor- 
tant especially when the orbits must be searched numer- 
ically. This is achieved by constructing cross- corr elated 
periodic orbit sums as introduced in Section |2.3| which 



are then harmonically inverted with the generalized filter- 
diagonalization method of Section |3.2[ The efficiency of 
the method will be discussed for various sets of opera- 
tors and various sizes of the cross-correlation matrix. It is 
also possible to include higher order h corrections in the 
cross-correlation signal which will allow us to calculate h 
corrections even for nearly degenerate states. 

Some concluding remarks are given in Section ra. 



2 Periodic orbit theory for integrable systems 

2.1 EBK quantization and Berry- Tabor formula 

Integrable systems are characterized by the property that 
their dynamics can be expressed in action-angle variables. 
The action variables, which are defined on certain "irre- 
ducible" paths, are constants of motion. In the 2n-dimen- 
sional phase space, the motion of an integrable system is 
restricted to n-dimensional tori, which are given by the 
values of the action variables. 

A well-established method for the semiclassical quanti- 
zation of integrable systems is the EBK torus quantization 
scheme, which was developed by Einstein, Brillouin and 
Keller ||,|,|. In the EBK theory, the energy eigenvalues 
of the system are directly associated with certain classical 
tori. These tori are defined by the EBK conditions, which 
select special sets from all possible values of the action 
variables of the system. Each such set corresponds to a 
quantum mechanical eigenstate of the system. The tori 
selected by the EBK conditions are usually not rational, 
i.e., the orbits on these tori are usually not periodic. 

For many physical systems the application of the EBK 
quantization scheme is a nontrivial task. Especially for 
non-separable or near-integrable systems the irreducible 
paths are difficult to find. Most importantly, as already 
discussed by Einstein Q| the torus quantization scheme 
cannot be extended to chaotic systems. For chaotic sys- 
tems, Gutzwiller derived a semiclassical expression for the 
density of states in terms of the periodic orbits of the cor- 
responding classical system: The semiclassical density of 
states consists of a smooth background and an oscillating 
part 

p(E)=p (E)+p° sc (E) (I) 



with 



p osc (E) = 



ST 



T 

± po 



7rft^r|det(M po -l)|Va 



( 2 ) 

The sum runs over all periodic orbits (po) of the system, 
including multiple traversals. Here, T and S are the period 
and the action of the orbit, M and p. are the Monodromy 
matrix and the Maslov index, and the repetition number 
r counts the traversals of the underlying primitive orbit. 

For integrable systems an analogous formula for the 
density of states in terms of a smooth part and oscillat- 
ing periodic orbit contributions was derived by Berry and 
Tabor ||. While in chaotic systems the periodic orbits 
are isolated, the periodic orbits of integrable systems are 
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all those orbits lying on rational tori - i.e., tori on which 
the frequencies of the motion are commensurable - and 
thus are non-isolated. The Berry- Tabor formula gives the 
density of states in terms of the rational tori: 



with 



p osc (E) 



p(E)=p (E)+p° sc (E) (3) 



cos^m/S- — \ira ■ M + t7t/3m) 



IMlK"- 1 ) |wm| |A-(I M )|* 



( 4 ) 

The sum runs over all rational tori at energy E, charac- 
terized by the frequency ratios given by the ray of integer 
numbers M. The sum includes cases where the Mj are 
not relatively prime, M = r/n, which corresponds to mul- 
tiple traversals of the primitive periodic orbits on the torus 
characterized by fi. Here, n is the dimension of the system, 
Im and o;m are the values of the action variables and the 
frequencies on the torus, 5m is the action of the periodic 
orbits on the torus, and K is the scalar curvature of the 
energy contour. The components of a are the Maslov in- 
dices of the irreducible paths on which the action variables 
are defined, and the phase (5 is obtained from the second 
derivative matrix of the action variables in terms of the 
coordinates. In contrast to the EBK torus quantization, 
there is no direct relation between the eigenvalues of the 
system and the tori which enter the Berry- Tabor formula. 

Both the EBK torus quantization and the Berry- Tabor 
formula are semiclassical theories delivering lowest order 
h approximations to the exact quantum eigenvalues. In 
general, the results of the two approaches can only be 
expected to be the same in lowest order of h but not nec- 
essarily beyond. However, it was shown in Ref. that 
for the circle billiard, which will be discussed in the fol- 
lowing sections, the two approaches are in fact equivalent 
and should yield exactly the same results. 

Eq. (Q) can be simplified for the special case n = 2, 
i.e., for two-dimensional regular systems |25 : 



n osc (F \ _ 1 \ " T M 

9 [ ] ^ l2 htMT\ g ^ 



M 



( 5 ) 

The sum runs over all rational tori at energy E, character- 
ized by the frequency ratio given by the integer numbers 
M = [M\,M-2), including multiple traversals (i.e., cases 
where M\,M% are not relatively prime). Here, Tm is the 
traversal time, and gE is the function describing the en- 
ergy surface: H{I\,l2 = gsih)) = E, where I\ and I 2 
are the action variables. The Maslov index r^M is obtained 
from the Maslov indices a\, a 2 of the paths on which the 
action variables are defined: 



Vm = (Afiai + M 2 a 2 ) - 9{g%] 



(6) 



where is the Heaviside step function. 

The semiclassical density of states can be expressed in 
terms of the semiclassical response function g(E): 



P(E) 



1 



Im g(E) . 



(7) 



For both chaotic and regular systems the response func- 
tion is of the form 



g{E) =g {E)+Y. A ^ Sva 



(8) 



po 



where the amplitudes are given by the Gutzwiller or the 
Berry- Tabor formula, respectively. 

In practical applications both the Gutzwiller formula 
(^|) and the Berry- Tabor formula (Q) suffer from the prop- 
erty that the periodic orbit sums usually do not converge. 
Depending on the system in question, this problem may 
be overcome, e.g., by convolution of the periodic orbit sum 
with a suitable averaging function |24j . But even then the 
convergence will usually be slow, and a large number of 
orbits has to be included in order to obtain the semiclas- 
sical eigenvalues. In the following sections, we will demon- 
strate how the convergence problem can be circumvented 
by the harmonic inversion method and the eigenvalues can 
be calculated from a relatively small number of periodic 
orbits. 



2.2 Application to the circle billiard 

We now apply the Berry- Tabor formula to the circle bil- 
liard as a specific separable system with two degrees of 
freedom. 

For completeness and comparisons with the results 
from periodic orbit theory we first briefly review the quan- 
tum mechanical expressions and the EBK quantization 
condition. The exact quantum mechanical energy eigen- 
values of the circle billiard with radius R are given by the 
condition 



J m {kR) = , m E Z 



E = 



h 2 k 2 

2M 



(9) 



where J m are the Bessel functions of integer order. Here, 
M den otes th e mass of the particle, E is the energy, and 
k = \/2ME/h is the wavenumber. The corresponding 
wave functions are given by 



ip(r,tp) = J m (kr)e l 



(10) 



As J_ m (x) = (—l) m J m (x), all energy eigenvalues be- 
longing to nonzero angular momentum quantum number 
(to 0) are twofold degenerate. In the following the ex- 
act quantum mechanical results for the circle billiard are 
used as a benchmark for the development and applica- 
tion of semiclassical quantization methods for integrable 
systems. 

The circle billiard problem in two dimensions is sepa- 
rable in polar coordinates. The semiclassical expressions 
for both EBK torus quantization and the Berry- Tabor for- 
mula for the density of states are based on the action-angle 
variables associated with the angular </?-motion and the 
radial r- motion (24lEq]. The action variables are given by 



(11) 



I v = — (f> pip dip = L 
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Ir = 



2n 



p r dr 



- ( ^2MER 2 -L 2 -\L\ 



arccos ■ 



\L\\ 



V2MERJ 



(12) 



where E and L are the energy and the angular momentum, 
respectively. Quantization of the action variables 



Iu, = [m-\ — r-) h , me 



H, rc = 0,l,2,. 



(13) 
(14) 



with OLtp = and a r ~ 3 for the circle billiard yields the 
EBK quantization condition 



v/(fci?) 2 - m 2 



Iml / 3\ 
\m\ arccos— = ln+ -) tt , (15) 



where L = mh are the angular momentum eigenvalues. 

The frequencies of the classical motion on the two- 
dimensional tori are given by 



dE 
dE_ 



2E 



^2MER 2 - L 2 

2nE 
V2MER 2 - L 2 



V2MER 



(16) 



(17) 



The Berry- Tabor formula includes all tori with a rational 
frequency ratio, i.e., tori on which the orbits are periodic. 
In the case of the circle billiard, the rational tori are given 
by the condition 

i 1 

(18) 



Mp 

M r 



with positive integers M r , M v and the restriction 



M r > 2M ir 



(19) 



The periodic orbits of the circle billiard have the form of 
regular polygons. The numbers M r and M v can be shown 
to be identical with the number of sides of the correspond- 
ing polygon and its number of turns around the center of 
the circle, respectively p7[ . Some examples are given in 
Figure [j]. A pair of numbers {M r ,M v ) which arc not 
relatively prime corresponds to multiple traversals of a 
primitive periodic orbit. 

The classical action of the periodic orbits is given by 



S M = 27rM„/,l M) + 27rM r 4 



M) 



' M<r 



V2MER 2M r sin | — p -k 

Mr 



As in all billiard systems, the action scales like 

S/h = ws , 
here with the scaling parameter 



w 



V2ME R/h=kR 



(20) 



(21) 



(22) 





(2,1) 



(3,1) 



(4,1) 





(5,2) (7,3) 

Fig. 1. Some examples of periodic orbits of the circle billiard. 
The orbits are labeled with the numbers (M r , M v ) which corre- 
spond to the number of sides of the polygons and the number of 
turns around the center. The angle 7 is given by 7 = nM v /M r . 



and the scaled action 



2M r sin 



Mp_ 

Mr 



(23) 



The form of the corresponding classical trajectory is in- 
dependent of w. For the circle billiard with unit radius 
R = 1, the scaling parameter w is identical with the 
wavenumber fc, and the scaled action is the length of the 
orbit. 

For the semiclassical density of states, we start from 
the special version of the Berry- Tabor formula presented 
in Eq. (j|). Using the relation 



PH = - P (E) 



(24) 



valid for billiard systems, we introduce the density of 
states depending on the scaling parameter w. Evaluating 
the different expressions in (|5|) for the circle billiard then 
finally leads to 



p osc (w) = Im g osc {w) 

7T 



with 



osc / \ 

9 (w) 



I — 3/2 
v M r 



,i(u;s M -|j\/ r 7r-f ) 



(25) 



(26) 



where we have used the relations a v = and a r = 3 for 
the Maslov indices. The sum runs over all pairs of positive 
integers M = (M r , M v ) with M r > 2M V . The degeneracy 
factor 



m M 



M r 

Mr 



= 2M {f 
> 2M K 



(27) 



accounts for the fact that all trajectories with M r ^ 2M V 
can be traversed in two directions. 
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(8,1) 



(12,1) 





(13,2) 



(17,2) 



(16,1) 



O 



(21,2) 



Fig. 2. Behavior of the periodic orbits of the circle billiard for 
large M r . The orbits are labeled with the numbers (M r ,M v ). 
As M r — » co at constant M v , the shape of the orbit converges 
against M v times the circular border of the billiard. The length 
of each single side of the orbits tends to zero, while the to- 
tal length of the orbit 2M r sm(jiM v /M r )R converges against 
2nM v R (R: radius of the billiard). The circle billiard there- 
fore possess accumulation points of orbits at scaled actions of 
multiples of 2n. 



Due to the rapid increase of the number of periodic or- 
bits with growing action, the sum ( p6| ) does not converge. 
In our case, the problem is even more complicated by the 
fact that there exist accumulation points of periodic or- 
bits at scaled actions of multiples of 2ir (see Fig. |J), which 
means that we cannot even include all periodic orbits up 
to a given finite action. In Ref. J24J] the convergence 
problem of the sum (^6|) was solved by averaging it with 
a Gaussian function. The semiclassical eigenvalues of the 
circle billiard were calculated from the periodic orbit sum 
by including a very large number of periodic orbits. Our 
aim is to demonstrate that by using the harmonic inver- 
sion scheme, we can obtain eigenvalues of w = kR from a 
relatively small number of periodic orbits. We will return 
to this problem in Section ^. 



2.3 Semiclassical matrix elements 

The semiclassical trace formula for both regular and 
chaotic systems can be extended to include diagonal ma- 
trix elements. The calculation of individual semiclassical 
matrix elements is an objective in its own right. Further- 
more, the extended trace formulae will allow us to con- 
struct cross-correlated periodic orbit signals and thus to 
significantly reduce the required number of orbits for peri- 
odic orbit quantization, as we will demonstrate in Section 



5.4 



Here, we will briefly recapitulate the basic ideas and 
equations. 

Both Gutzwiller's and Berry and Tabor's formula give 
the semiclassical response function as a sum over contri- 
butions from periodic orbits (see Eq. (||)). The quantum 
mechanical response function is the trace over the Green 



function G 



E- 



iO 



tr G+ 



(28) 



As a generalization, one can consider the quantum me- 
chanical response function weighted with the diagonal ma- 
trix elements of some operator A, i.e., 



9 A,. 



^) = E 



(n\A\n) 



E — E„ 



iO 



= tr (G+A) 



(29) 



The semiclassical approximation to the extended response 
function (29) is obtained by weighting the contributions 
of the periodic orbits in (g) with the average A p of the 
corresponding classical quantity A(q, p) over the periodic 
orbits: 



g A (E) = g A ,o(E) + ^ A po A p e* ' 



(30) 



For chaotic systems, the average is taken over one period 
T p of the isolated periodic orbit ^,^9,30 : 



A, 



T, 



A(q(t),p(t)) dt 



(31) 



p Jo 



For an iV-dimensional integrable system, the quantity A 
has to be expressed in action-angle variables (I, 6) and 
averaged over the rational torus pq| : 



At, 



1 



(27T) 



A{i,e)d N e . 



(32) 



Eq. ( p9| ) can even be further generalized by introducing 
a second operator B and considering the quantity [|lq] 



(E) 



E 



(n\A\n) (n\B\n) 



E — E„ 



iO 



(33) 



If either Aov B commutes with the Hamiltonian, Eq. ( |3^ ) 
can be written as a trace formula and a calculation similar 
to that in Efl] yields the semiclassical approximation 



9ab(E) = g A B,o( E ) + E ApoA p B p e 



(34) 



po 



Note that for general operators A and B, Eq. ( |33| ) cannot 
be written as a trace any more. However, strong numer- 
ical evidence was provided (for both regular and chaotic 
systems) that Eq. ([m]) is correct in general, i.e., even if 
neither operator A nor B commutes with the Hamiltonian 
[ ^6[ . For chaotic systems, a mathematical proof of Eq. (Q) 
is given in Ref. |18| . An analogous rigorous derivation for 
integrable systems is, to our knowledge, still lacking. 

In Refs. |T^,|l8|, the relations ( p3| ) and (34) were gen- 
eralized to products of diagonal matrix elements of more 
than two operators. As a further extension, we can also 



6 Kirsten Weibert et al.: Use of harmonic inversion techniques in the periodic orbit quantization of integrable systems 



introduce functions of diagonal matrix elements in the re- 
sponse function: 



9f(A),qm(E) 



V f((n\A\n)) 



(35) 



By a Taylor expansion of the (sufficientlysmooth) func- 
tion / and using the results of Refs. [ p^|]T^ | for multiple 
products of matrix elements, we obtain the semiclassical 
approximation 



9f(A)(E) =<?/ ( a),o(£) +J2 A pof(A p )e 



(36) 



po 



We will use the extended trace formulae in combination 
with an extension of the harmonic inversion procedure to 
cross-correlated signals in order to significantly reduce the 
number of orbits which have to be included in the periodic 
orbit sum. 

The diagonal matrix elements obtained from the ex- 
tended trace formulae are semiclassical approximations to 
the exact quantum matrix elements. For the circle bil- 
liard, we can compare these values to those given by EBK 
theory. According to EBK theory, the diagonal matrix ele- 
ment of an operator A with respect to an eigenstate |n) is 
obtained by averaging the corresponding classical quantity 
A(l, 6) over the quantized torus related to this eigenstate: 



(n\A\n) 



1 



(2tt 



A{l ni 6 n )d N 6, 



(37) 



Note the difference to Eq. (|32|), where the average is taken 
over the rational tori. 



developed by Gaspard and Alonso [|31|,|32|] and by Vat- 
tay and Rosenqvist |3^,^J,^5|, following two different ap- 
proaches. Vattay and Rosenqvist compute the corrections 
by solving the local Schrodinger equation in the neighbor- 
hood of periodic orbits. They introduce a quantum gen- 
eralization of the Gutzwillcr formula which contains these 
local eigenvalues. The results of Refs. ||^,[54],[35| cannot di- 
rectly be applied to integrable systems, as the derivations 
are valid only for isolated periodic orbits. To our knowl- 
edge, a general theory for h corrections to the Berry- Tabor 
formula does not yet exist. 

Nevertheless, for the circle billiard we have succeeded 
in obtaining an explicit expression for the first order h 
corrections to the Berry- Tabor formula. The calculations 
are quite lengthy and are therefore deferred to Appendix 
|A|. Our final result for the first order H amplitude of the 
circle billiard in (^9|) reads: 



A {1) 



1 



3 sin 7 6 sin 7 



(40) 



with 7 = nM v /M r and Apo the zeroth order amplitudes 
given by the Berry- Tabor formula. Using the zeroth order 
amplitudes from Eq. (Eq), we finally obtain 



m , — — -2 sin 2 7 — 5 
6sm J/ ^7 



i(|M r 7r-f) 



(41) 



As explained above our derivation of Eq. (|4l|) cannot be 
applied to general integrable systems. It will be an inter- 
esting task for the future to develop a general theory for 
the higher order h corrections to the Berry- Tabor formula. 



2.4 Higher order h corrections 

The Berry- Tabor formula for integrable systems and 
Gutzwiller's trace formula for chaotic systems are only the 
leading order terms of an expansion of the density of states 
in terms of H. In billiard systems, the scaling parameter 
w of the classical action (cf. Eq. ([H])) is proportional to 
ftr 1 and thus plays the role of an inverse effective Planck 
constant, 

w = n;i . (38) 

The h expansion of the response function can therefore be 
written as a power series in terms of w _1 pa]: 



3 Harmonic inversion by filter-diagonalization 

The quantization of the periodic orbit sum as well as the 
analysis of quantum spectra in terms of the periodic or- 
bits can be reformulated as a harmonic inversion problem 
of formulae which have been introduced in the previous 
Sections. Before discussing these applications in Sections 
[| and [)|~we will now briefly recapitulate the basic ideas 
and the technical tools of har monic inversion by filter- 
diagonalization. In Section 3.1 we will start with the har- 
monic inversion of a sin gle f unction. The equations will 
be generalized in Section 3.2 to the harmonic inversion of 
cross-correlated signals. 



n=0 



n=0 



The zeroth order amplitudes A v °o are those of the Berry- 
Tabor or Gutzwillcr formula, respectively, whereas for 
n > 0, the amplitudes A$ give the n th order corrections 
g n (w) to the response function. Explicit expressions for 
the first order correction terms for chaotic systems were 



(39) 



3.1 Harmonic inversion of a single function 

The harmonic inversion problem can be formulated as a 
nonlinear fit of a signal C(t) to the form 



C(t) 



k 



fee 



(42) 
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where dk and u>k are generally complex variational param- 
eters. Other than, e.g., in a simple Fourier transformation 
of the signal, there is no restriction to the closeness of 
the frequencies ujk ■ Solving ([l2]) will therefore yield a high 
resolution analysis of the signal C(t). The signal length 
required for resolving the frequencies uik by harmonic in- 
version can be estimated to be 



The frequencies u>k are then obtained by solving the gen- 
eralized eigenvalue problem 



U (p) B fe = e - ip ™ fc Bt 



(50) 



The amplitudes dk can be calculated from the eigenvectors 
and are given by 



47rp(w), 



(43) 



M 



where p(ui) is the mean density of frequencies in the range 
of interest. 

A method which has proven very useful for solving the 
harmonic inversion problem is the filter-diagonalization 
procedure |H^,^(|. This procedure allows us to compute 
the frequencies Uk in any small interval [w m i n , w max ] given. 
The idea is to consider the signal C(t) on an equidistant 
grid 

c n = C(nr); n = 0,1,2,... (44) 

and to associate c n with an autocorrelation function of a 
suitable fictitious dynamical system, described by a com- 
plex symmetric effective Hamiltonian H c r: 



rH, 



(45) 



Here, the brackets denote a complex symmetric inner 
product (a\b) = (b\a), i.e., no complex conjugation of ei- 
ther a or b. The harmonic inversion problem can then be 
reformulated as solving the eigenvalue problem for the ef- 
fective Hamiltonian H c g . The frequencies Uk are the eigen- 
values of the Hamiltonian 



H eS \Tk) = uJk\T k ) , 



(46) 



and the amplitudes are obtained from the eigenvectors T^: 



dk = (#o|Tfc) 



(47) 



The filter-diagonalization method solves this eigenvalue 
problem in a small set of basis vectors S^ . The Hamil- 
tonian and the initial state <Pq do not have to be known 
explicitly but are given implicitly by the quantities c„. In 
detail, the procedure works as follows: 

A small set of values tpj in the frequency interval of 
interest is chosen. The set must be larger than the number 
of frequencies in this interval. The values ipj are used to 
construct the small Fourier-type basis 



A I 



(48) 



n=0 



The matrix elements of the evolution operator at a given 
time pr in this basis can be expressed in terms of the 
quantities c„: 



M M 



n=0 n'=0 



(49) 



dk = B jk ^2 



(51) 



The values of ujk and dk obtained by the above procedure 
should be independent of p. This condition can be used 
to identify non-converged frequencies by comparing the 
results for different values of p. The difference between 
the frequency values obtained for different p can be used 
as a simple error estimate. 



3.2 Harmonic inversion of cross-correlated signals 

An important extension of the filter-diagonalization 
method for harmonic inversion is t he g eneralization to 
cross-correlation functions [p^[^ ,^,|l7|\ This extended 
method allows us to significantly reduce the signal length 
required to resolve the frequencies contained in the signal. 
The idea is not to consider a single signal C(t) as given in 
Eq. ( f42|) but a set of cross-correlated signals 



Caa'ip) ^ y d aa > ^k& 



with the restriction 



—iuikt 



a,a' = l,...,N (52) 



d aa ',k = b a ,kb a ',k ■ (53) 

This set of signals considered on an equidistant grid 

Cnaa' = C aa >{nT) \ 71 = 0,1,2,... (54) 

is now associated with a time cross-correlation function 
between an initial state ^> a and a final state <S> a i : 



(55) 



Again, the frequencies LUk are obtained as the eigenvalues 
of the effective Hamiltonian H c s . The amplitudes are now 
given by the relation 



i$a\Tk) 



(56) 



In analogy to the procedure described in Section 3.1, this 
eigenvalue problem is solved in a small set of basis vectors 
tf'ja in order to obtain the frequencies in a given interval 

The advantage of the above procedure becomes evi- 
dent if one considers the information content of the set 
of signals. Due to the restriction (^3|), the N x N set of 
signals C aQ < (t) may contain N independent signals, which 
are known to possess the same frequencies u>k- This means 



8 Kirsten Weibert et al.: Use of harmonic inversion techniques in the periodic orbit quantization of integrable systems 



that, at constant signal length, the matrix may contain 
N times as much information about the frequencies as 
a single signal, provided that the whole set is inverted si- 
multaneously. On the other hand, the information content 
is proportional to the signal length. This means that the 
signal length required to resolve the frequencies in a given 
interval is reduced by a factor of N. This statement clearly 
holds only approximately and for small matrix dimensions 
N. However, a significant reduction of the required signal 
length can be achieved. 



4 High resolution analysis of quantum spectra 

Harmonic inversion is a powerful tool to calculate the clas- 
sical periodic orbit contributions to the density of states 
from the quantum mechanical eigenvalues or from experi- 
mental spectra, thus delivering a high resolution analysis 
of the spectra in terms of the classical orbits. The method 
allows us, e.g., to resolve clusters of orbits or to discover 
hidden ghost orbit contributions in the spectra, which 
would not be resolved by conventional Fourier analysis 
of the spectra Jl3|,p3|. Here, we will analyze the quantum 
spectra of the circle billiard as a representative of inte- 
grable systems. The analysis will verify the validity of the 
Berry- Tabor formula and its extensions to semiclassical 
matrix elements and higher order h corrections discussed 
in Section 0. 



4.1 Leading order periodic orbit contributions to the 
trace formula 



4.1.1 General procedure 

In this section we develop the general procedure for the 
analysis of quantum spectra in terms of periodic orbits 
by harmonic inversion. This procedure is universal in the 
sense that it can be applied to both regular and chaotic 
systems. We will apply it to the circle billiard as a repre- 
sentative of regular systems in the next section. 

We start from the semiclassical density of states given 
by the Berry- Tabor or the Gutzwiller formula. As in Sec- 
tion ||, we consider scaling systems where the density of 
states depends on the scaling parameter w [w — kR for 
the circle billiard], i.e., p(w) = — (l/7r)Im g(w) with g(w) 
the semiclassical response function. Both the Berry- Tabor 
and the Gutzwiller formula give the oscillating part of the 
response function in the form 

po 

where the sum runs over all rational tori (regular systems) 
or all periodic orbits (chaotic systems) of the underlying 
classical system, respectively. Here, S po is the action of the 



periodic orbit. The form of the amplitude A po depends on 
whether the system is chaotic or regular and also contains 
phase information. 

The exact quantum mechanical density of states is 
given by 

Pqm{w) = ^TOfc(5(w - Wk) , (58) 
k 

where the Wk are the exact quantum eigenvalues of the 
scaling parameter and the rrik are their multiplicities. The 
analysis of the quantum spectrum in terms of periodic 
orbit contributions can now be reformulated as adjusting 
the exact quantum mechanical density of states ( p)8| ) to 
the semiclassical form 

P osc (w) = - ilm g osc (w) 

= ~E {A po e iws »° - A; o e- l ^°) . (59) 

po 

If the amplitudes A po do not depend on w, the semiclas- 
sical density of states is of the form (^) [here, with w 
playing the role of t and s po playing the role of u>k\- That 
means, we have reformulated the problem of extracting 
the periodic contributions from the quantum spectrum as 
a harmonic inversion problem. In the fitting procedure, 
we ignore the non-oscillating, smooth part of the density 
of states. This part does not fulfill the ansatz (E5|) of the 
harmonic inversion method and therefore acts as a kind 
of noise, which will be separated from the oscillating part 
of the "signal" by the harmonic inversion procedure. 

In practice, in order to regularize the S functions in 
(p8[), we convolute both expressions (|5^) and ( |59| ) with a 
Gaussian function, 

C a {w) = / p(w')e- {w - w '^^ a2 dw' . (60) 

V27TCT J_oo 

In our calculations, we usually took the convolution width 
to be about three times the step width r in the signal (fli]) . 
Typical values are r = Aw = 0.002 and a = 0.006. The 
resulting quantum mechanical signal is 

C qm , CT M - -J=- V m k e-^- w ^' 2 ° 2 , (61) 
V 2na 

and the corresponding semiclassical quantity reads 
C a {w) =-iV (A po e tws ?° - A*e~* ws *°) e~^ s l° . 

po 

The above procedure still works if the amplitudes in 
( p7| ) are not independent of w but possess a dependency 
of the form 

A po = w a a po , (63) 

which is, for example, the case for regular billiards. This 
dependency can be eliminated |l5|] by replacing the semi- 
classical response function g(w) with the quantity 

g '( w ) = w- a g(w) = w- a g Q (w) + £ a po e lws *° . (64) 

po 
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When introducing the corresponding quantum mechanical 
response function 



3qmM = ' 



m k w k 



Wk 



iO 



(65) 



the procedure can be carried out for p (w) = 
(—l/ir)Img'(w) as described above. 

In addition to cons ider ing the pure density of states, 
the relations of Section 2.2 can be used to calculate the av- 



erages of classical quantities over the periodic orbits from 
the quantum diagonal matrix elements of the correspond- 
ing operators. If we start from the extended quantum re- 
sponse function (29), including diagonal matrix elements 
of some operator A, the analysis of the signal should again 
yield the actions s po as frequencies but with the ampli- 
tudes weighted with the classical averages A p of the corre- 
sponding classical quantities. In the same way, we can also 
use the extended response function (|3^), which includes 
diagonal matrix elements of two different operators. 



4.1.2 Application to the circle billiard 

For the circle billiard, the oscillating part g osc (w) of the 
semiclassical response function is given by Eq. (|26|). If 
one eliminates the dependency of the amplitudes on w 
by defining 

p'H = -^pH , (66) 



gi(|M r 7r-f ) e -iwstA 



the resulting expression for the density of states 

o 3/2 

> mivr - 

1V1 

e -i( | M r7 r- f ) e iwsw ] (57) 



is of the form (|42|), here with Sm. playing the role of Wfc. 
The quantum mechanical quantity corresponding to ( |66[ ) 

TO 

P'qm( W ) = -?=~ S ( W ~ W k) ■ (68) 



IS 



In addition to analyzing the pure quantum spectrum 
of the circle billiard, we also considered spectra weighted 
with di agon al matrix elements of different operators (cf. 
Section 2.3). We used three different operators, viz. 



the absolute value of the angular momentum L as an 
example of a constant of motion, 
the distance r from the center as an example of a quan- 
tity which is no constant of motion, 



the variance of the radius (r 



as an example 



using the relation ( p4| ) for products of operators. 

The classical angular momentum L is proportional to w, 
which means that when constructing the signal for L, g(w) 
now has to be multiplied by u>~ 3 / 2 instead of u> -1 / 2 (cf. 
Eq. ©). 
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— 
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Ilk 
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16 
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20 



Fig. 3. Periodic orbit contributions to the trace formula calcu- 
lated from the quantum spectrum of the circle billiard, includ- 
ing different operators. The quantities plotted are the classical 
amplitudes Am. = m-MS M 2 /M% times the classical averages of 
the operators indicated, versus the scaled actions of the or- 
bits. Crosses: values obtained from the quantum spectrum by 
harmonic inversion. Solid lines: values obtained from classical 
calculations. L: angular momentum in units of hw (w: scaling 
parameter), r: distance from the center in units of the radius 
R of the billiard. 



We calculated the scaled actions and classical ampli- 
tudes of the periodic orbits in the interval sm € [15,23]. 
The signal was constructed from the exact zeros of the 
Bessel functions, up to a value of w max — 500. The ac- 
curacy of results is improved if we cut off the lower part 
of the signal, using only zeros larger than w m i n = 300. A 
possible explanation for this is that the low zeros are in a 
sense "too much quantum" for the semiclassical periodic 
orbit sum. 

Figure || shows the results of our calculation. The po- 
sitions of the solid lines are the scaled actions of the clas- 
sical periodic orbits, their heights are the classical ampli- 
tudes mMS^ 2 /M 2 times the respective averaged classical 
quantity. The crosses are the results obtained by harmonic 
inversion of the signal constructed from the zeros of the 
Bessel functions. There is an excellent agreement between 
the spectra, illustrating the validity of the Berry- Tabor 
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formula and its exten sion to semiclassical matrix elements 
discussed in Section 2.3, The examined interval contains 
an accumulation point of orbits (s = 6tt). Here, only those 
orbits were resolved which were still sufficiently isolated. 

It might be surprising that, although the Berry- Tabor 
formula only gives a semiclassical approximation to the 
density of states and we started from the exact quantum 
mechanical density, our results for the periodic orbit con- 
tributions are exact and do not show any deviations due to 
the error of the semiclassical approximation. The reason 
for this will become obvious in the following Section. 




200 200.1 200.2 200.3 200.4 200.5 200.6 



4.2 Higher order h corrections to the trace formula 

An interesting a ppli cation of the method described in the 
previous Section 4.1 is the determination of higher order h 



contributions to the periodic orbit sum. The higher orders 
can be obtained by analysis of the difference spectrum 
between the exact quantum and semiclassical eigenvalues, 
as we will show below. 



As explained in Section 2.4, the Berry- Tabor formula 
for integrable systems as well as the Gutzwiller formula 
for chaotic systems are the leading order terms of an ex- 
pansion of the density of states in terms of h. For scaling 
systems, this expansion can be put in the form (cf. 



n=0 



E 

n=Q 



1 



po 

(n) 



(69) 



i 



Provided that the amplitudes A po in fl69| ) do not de- 
pend on w, only the zeroth order term fulfills the ansatz 
( fiHI ) for the harmonic inversion procedure with constant 
amplitudes and frequencies. In systems like regular bil- 
liards, where the amplitudes possess a w dependence of 
the form — w a a p r $ , the same argumentation holds if 
we co nsider g'(w) = w~ a g(w) instead of g(w) (cf. Section 
l.lj ) . This is the reason why the analysis of the quantum 

spectrum yields exactly the lowest order amplitudes A$ , 
without any deviations due to the semiclassical error: As 
the higher order terms do not fulfill the ansatz, the Apo 
are the best fit for the amplitudes. The higher oder terms 
have similar properties as a weak noise and are separated 
from the "true" signal by the harmonic inversion proce- 
dure. 

Although the direct analysis of the quantum spectrum 
only yields the lowest order amplitudes, higher order cor- 
rections can still be extracted from the spectrum by har- 
monic inversion. Assume that the exact eigenvalues Wk 
and their (n — l) st order approximations Wk, n -i are given. 
We can then calculate the difference between the exact 
quantum mechanical and the (n— l) st order response func- 
tion 



qm 



n—l oc oo 1 

h - E ftH - E »(») = E ^ £4£ 



]=n 



j=n po 



(70) 



Fig. 4. Part of the difference spectrum Ap qm (w) — Apebk(w) 
between the exact quantum mechanical and the semiclassical 
density of states. The absolute values of the peak heights mark 
the multiplicities of the states. 



The leading order terms in 
cation with w n yields 



Oh are 



w 



multipli- 



10 



s qn » 



n-l 

E 

3=0 



9j( w ) 



E4"V 



o 



-n na (71) 

In (J7IJ) we have restored the functional form ([I2j). The 

harmonic inversion of the function ([71]) will now provide 



the periods Spo and the n th order amplitudes Ap^ of the 
h expansion (p9|). 

In p ractic e, we will follow the procedure outlined in 
Section [iLl] to construct a smooth signal, i.e., we consider 
the densities of states p(w) = — (l/7r) Jjoag(w) rather than 
the response functions g(w), and smoothen the signal by 
convoluting it with a Gaussian function. 

For the circle billiard, the exact quantum eigenval- 
ues are given by the condition (||), while the zeroth or- 
der eig envalues are equ al to the EBK eigenvalues given 
by (|1^) (cf. Section 2A ) . From the difference between the 
exact and the semiclassical density of states, we can cal- 
culate the amplitudes Ap 1 } of the first order correction to 
the trace formula. 

We analyzed the difference spectrum between exact 
and EBK eigenvalues of the circle billiard in the range 
100 < w < 500. Figure [| shows a small part of this dif- 
ference spectrum. The results of the harmonic inversion 
of the spectrum are presented in Figure | The crosses 
mark the values 



/(7) 



y/wM r 



I — \ Apo I ' 



(72) 



with 7 = TtM v /M r which we obtained for the periodic 
orbits by harmonic inversion of the difference spectrum. 
The crosses are labeled with the numbers (M r , M v ) of the 
orbits. The solid line in Fig. || is the theoretical curve 



/(7) 



2 sin 7 



3 sin' 



3/2 



7 



(73) 
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Fig. 5. First order h corrections to the Berry- Tabor formula 
for the circle billiard. Crosses: amplitudes (2/ '^/nAlr) \ A^p] \ 
obtained by harmonic inversion. The orbits are labeled by 
the numbers (M r ,M<p). Solid line: theoretical curve f("f) = 
(5 -2 sin 2 7)/(3sin 3/2 7 ). 



This fitting problem cannot be solved directly, but can be 
reformulated as a harmonic inversion problem jll],[l^] . The 
first step of the reformulation is a Fourier transformation 
of the response functions with respect to w: 



1 f°° 
C(s) = — / g(w)t 



,J dw 



(76) 



In the semiclassical response function, we only consider 
the oscillating part of g(w). The smooth part, which does 
not possess a suitable form for the harmonic inversion 
method, would only give a contribution to the signal for 
very small s. Assuming that the amplitudes in ( |74| ) do not 
depend on w, the result of the Fourier transformation is 



C(s) = ^A po S(s - s po ) 

po 

C qm (s) = -iJ2 m ke~ isw 



(77) 
(78) 



which results from our analytical express ion ( J4l| ) for the 
first order amplitudes discussed in Section 2.4. The results 
obtained by harmonic inversion are in excellent agreement 
with the theoretical curve, which clearly illustrates the 
validity of Eq. (p|). 



5 Periodic orbit quantization 



5.1 General procedure 

We now turn to the problem of extracting eigenvalues from 
the periodic orbit sum. We will demonstrate that the har- 
monic inversion procedure, which has already been suc- 
cessfully applied to extract the eigenvalues of chaotic sys- 
tems | p[p^ , can be used for integrable systems as well 
when starting from the Berry- Tabor formula. 

As previously (see Section ^) , we consider scaling sys- 
tems and start from the response function 



go (w) + } ^A po e~ 



po 



(74) 



depending on the scaling parameter w. The amplitudes 
„4 po are those of the Berry- Tabor or the Gutzwiller for- 
mula for regular and chaotic systems, respectively. The 
periodic orbit sum in ( [74]) usually does not converge, or, 
at least, the convergence will be very slow. In practice, 
especially for chaotic systems, only the orbits with small 
scaled actions will be known. Nevertheless, the eigenval- 
ues of the scaling parameter can still be extracted from the 
periodic orbit sum. The central idea is to adjust Eq. ([74]), 
with the sum including periodic orbits up to a finite action 
s max , to the functional form of the corresponding quantum 
mechanical response function 



w — Wk + i0 



(75) 



Like in Section |4.1.1| we convolute the signals ( f77| ) and 
(f7q) with a Gaussian function with width er, resulting in 



C a {s) 



-(S-Spo) 2 /^ 2 



po 

G,,rvx(«) = -*X! mfee r 
k 



(79) 
(80) 



Typical values of the convolution width are a = 0.006 for 
signals with step width As — 0.002. The eigenvalues of 
the scaling parameter are now obtained by adjusting the 
signal C a (s) to (|80|), which is of the functional form (p2|). 
The frequencies w k obtained by harmonic inversion ofthc 
signal (J73) are the eigenvalues of the scaling parameter 
w; from the amplitudes d k , the multiplicities can be 
calculated. 

Like the gener al pro cedure for analyzing quantum 
spectra (see Section 4.1.1 ) , the above procedure still works 
if the amplitudes in (|74|) are not independent of w but pos- 
sess a dependency of the form 



•/4,po — W dpo 



(81) 



We can again eliminate this dependency by replacing g(w) 
with the quantity 

g'(w) = w- a g(w) . (82) 
The semiclassical signal is then given by 



-( S -s po ) 2 /2<t 2 



(83) 



and the corresponding quantum mechanical signal reads 
C qm , ff (s) = -i J2 m k w^ a e-^ 2 <e^ . (84) 
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5.2 Semiclassical eigenvalues of the circle billiard 



5.2.1 Construction of the periodic orbit signal 

The semiclassical response function of the circle billiard is 
given by Eq. (Sfjl) . The amplitudes in (|26|) are proportional 



to w 1 / 2 . As described in Section 5.1, we eliminate this 



dependency on w by introducing the quantity 

g'{w) = w-^ 2 g(w) . (85) 



Applying Eqs. ( p3| ) and (84), we now obtain the semiclas- 
sical and the corresponding exact quantum signal for the 
circle billiard: 



2a 



3/2 

M 



e -iiM r 7T e -(s-s M y/2^ 



Cqm,cT (s) 



E 



2° w kp~ lsw k 



w. 



(86) 
(87) 



Eq. (|87|) possesses the functional form (B2) with 



, . m k _i 2 2 

d k = i _ < " *• 



Applying the harmonic inversion method to the signal ( p6[ ) 
should yield the eigenvalues of w as frequencies, with the 
amplitudes given by Eq. (|S8|). 



5.2.2 Results for the lowest eigenvalues 

We calculated the eigenvalues of the scaling parameter 
w = kR for the lowest states of the circle billiard from a 
signal of length s max = 150. For the construction of the 
signal, we chose a minimum length for the sides of the pe- 
riodic orbits as cut-off criterion at the accumulation points 
(cf. Fig. ||). We observed that the results were nearly in- 
dependent of the choice of this parameter, as long as the 
minimum length was not chosen too large. 

Table [l] presents the semiclassical eigenvalues Whi and 
multiplicities m^i obtained by harmonic inversion of the 
periodic orbit signal (|8^). For comparison, the exact quan- 
tum mechanical and the EBK results are also given in 
Table The eigenvalues obtained by harmonic inversion 
clearly reproduce the EBK eigenvalues within an accu- 
racy of 1CP 4 or better. The deviation of the ium from the 
EBK eigenvalues is significantly smaller than the error 
of the semiclassical approximation. The improvement of 
the semiclassical quantization by including higher order h 
corrections to the periodic orbit sum will be discussed in 



Section 5.3 



Note that for calculating the eigenvalues of the circle 
billiard by a direct evaluation of the periodic orbit sum, 
a huge number of periodic orbit terms is required, e.g., 
orbits with maximum length s max = 30 000 were included 



Table 1. Lowest eigenvalues w k and multiplicities of the 
scaling parameter w — kR of the circle billiard. w cx and ra ex : 
exact quantum values, webk: EBK eigenvalues, w^i and m^: 
values obtained by harmonic inversion of a signal of length 
s max = 150. The numbers n and m are the radial and angular 
momentum quantum numbers. The nearly degenerate eigen- 
values at « 11.0 and ~ 13.3 were not resolved. 



n 


m 


™ cx 


webk 






m hi 








2.404826 


2.356194 


1 


2.356204 


1.0005 





1 


3.831706 


3.794440 


2 


3.794444 


1.9983 





2 


5.135622 


5.100386 


2 


5.100391 


1.9996 


1 





5.520078 


5.497787 


1 


5.497785 


0.9988 





3 


6.380162 


6.345186 


2 


6.345191 


2.0000 


1 


1 


7.015587 


6.997002 


2 


6.997006 


2.0001 





4 


7.588342 


7.553060 


2 


7.553065 


1.9992 


1 


2 


8.417244 


8.400144 


2 


8.400149 


1.9998 


2 





8.653728 


8.639380 


1 


8.639404 


0.9987 





5 


8.771484 


8.735670 


2 


8.735677 


2.0013 


1 


3 


9.761023 


9.744628 


2 


9.744632 


1.9999 





6 


9.936110 


9.899671 


2 


9.899675 


1.9999 


2 


1 


10.173468 


10.160928 


2 


10.160932 


2.0000 


1 


4 


11.064709 


11.048664 


2 


11.048968 


4.0012 





7 


11.086370 


11.049268 


2 


2 


2 


11.619841 


11.608251 


2 


11.608256 


2.0006 


3 





11.791534 


11.780972 


1 


11.780978 


1.0001 





8 


12.225092 


12.187316 


2 


12.187319 


1.9993 


1 


5 


12.338604 


12.322723 


2 


12.322724 


2.0000 


2 


3 


13.015201 


13.004166 


2 


13.004168 


1.9997 


3 


1 


13.323692 


13.314197 


2 


13.315045 


4.0287 











9 


13.354300 


13.315852 


2 


1 


6 


13.589290 


13.573465 


2 


13.573465 


2.0000 


2 


4 


14.372537 


14.361846 


2 


14.361849 


1.9994 





10 


14.475501 


14.436391 


2 


14.436395 


2.0006 


3 


2 


14.795952 


14.787105 


2 


14.787076 


1.9909 


1 


7 


14.821269 


14.805435 


2 


14.805453 


2.0066 


4 





14.930918 


14.922565 


1 


14.922569 


1.0001 



in Ref. [Q. We obtained similar results using only or- 
bits up to length s max = 150. This demonstrates the high 
efficiency of the harmonic inversion method in extracting 
eigenvalues from the periodic orbit sum. The efficiency can 
even be further increased with the help of cross-correl ated 
periodic orbit sums as will be demonstrated in Section |5.4| . 

In Tabic [l] the exact multiplicities m cx of eigenvalues 
and the multiplicities mhi obtained by harmonic inversion 
also agree to very high precision. The deviations are one 
or two orders of magnitude larger than those in the fre- 
quencies, which reflects the fact that, with the harmonic 
inversion method using filter-diagonalization, the ampli- 
tudes usually converge more slowly than the frequencies. 

In the frequency interval shown, there are two cases 
of nearly degenerate frequencies which have not been re- 
solved by harmonic inversion of the periodic orbit signal 
with s max = 150. The harmonic inversion yielded only one 
frequency, which is the average of the two nearly degen- 
erate ones, with the amplitudes of the two added. These 
nearly degenerate states can be resolved when the signal 
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length is increased to about s max = 500 or with th e he lp 
of cross-correlated periodic orbit sums (see Section 5.4). 



5.2.3 Semiclassical matrix elements 

Using the extended periodic orbit sums discussed in Sec- 
tion p.3| , we can now also calculate semiclassical diagonal 
matrix elements for the circle billiard. 

Following the procedure described in Section 



5.1 



semiclassical signal can be constructed from the extended 
response function, the analysis of which should again yield 
the eigenvalues Wk as frequencies but with the amplitudes 
weighted with the diagonal matrix eleme nts. A s examples, 
we used the same operators as in Section 1.1.2 to calculate 
the diagonal matrix elements (L) , (r) , and the variance of 
the radius, Jr 2 ) — (r) 2 . 

Figure g shows the results in the range 25 < w < 30. 
For comparison, Fig. ^a presents the spectrum for the 
identity operator. The positions of the solid lines are the 
EBK eigenvalues, their heights are the semiclassical ma- 
trix elements obtained from EBK theory times the mul- 
tiplicities. The crosses are the results of the harmonic in- 
version of a signal of length s max = 300. The diagrams 
show excellent agreement between the results obtained by 
harmonic inversion and EBK torus quantization. This is 
even the case for the variance of r, which is a very small 
quantity. 

For the states shown in Fig. ^ we have also compared 
the semiclassical to the exact quantum mechanical matrix 
elements. The agreement is also excellent. The deviations 
between the semiclassical and quantum matrix elements 
are typically of the order of ~ 10~ 3 , which can well be 
understood by the semiclassical approximation. 



5.3 Higher order h corrections 

The eigenvalu es o f the circle billiard obtained in the pre- 
vious Section 5.2 are not the exact quantum mechanical 
eigenvalues but semiclassical approximations for the rea- 
son that the Berry- Tabor and the Gutzwiller formula are 
only the leading order terms of an exp ansi on of the den- 
sity of states in terms of H (see Section 2.4). We will now 
demonstrate how to obtain corrections to the semiclassi- 
cal eigenvalues from the h expansion (|39| ) of the periodic 
orbit sum 
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Fig. 6. Circle billiard: Results of the harmonic inversion of the 
semiclassical signal constructed from the periodic orbits sum 
including different operators. The quantities plotted are the 
diagonal matrix elements of the operators indicated times the 
multiplicities mk, versus the eigenvalues Wk of the scaling pa- 
rameter. Crosses: results of the harmonic inversion procedure. 
Solid lines: semiclassical matrix elements obtained from EBK 
theory and EBK eigenvalues. L: angular momentum in units 
of h, r: distance from the center in units of the radius R of the 
billiard. 



simplicity, we will assume in the following that the ampli- 
tudes A$ in ( |39| ) do not depend on w. Again, in systems 
where the amplitudes possess a w dependence of the form 



A (n) 



» 



the same line of arguments holds if we 
= w~ a g(w) instead of g(w) (cf. Section 



For periodic orbit quantization the zeroth order contri- 
butions A^ oo a re usually considered only. As demonstrated 



con sider g'(w) 

IB. 



in Section 5J (see Eqs. ([77]) and (|7§|) ) , the Fourier trans- 
form of the principal periodic orbit sum 



with w 
Eq. (|3 



h e fi an effective inverse Planck constant (see 



). The amplitudes A^?} are those of the Berry- 
Tabor or Gutzwiller formula. For n > 0, the amplitudes 
^Ipo'' (including also phase information) give the n th order 
corrections g n (w) to the response function g osc (w). For 



a ( s ) = ^4° o ^( s - Spo ) 

po 



is adjusted by application of the harmonic inversion tech- 
nique to the functional form of the exact quantum expres- 
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Cqm(s) 



E 



m k e 



with {wk,rrik} the eigenvalues and multiplicities. 

For n > 1, the asymptotic expansion ( |39| ) of the semi- 
classical response function suffers from the singularities 
at w = 0. It is therefore not appropriate to harmonically 
invert the Fourier transform of ( p9| ) as a whole, although 
the Fourier transform formally exists. This means that 
the method of periodic orbit quantization by harmonic 
inversion cannot straightforwardly be extended to the h 
expansion of the periodic orbit sum. Instead, we will cal- 
culate the correction terms to the semiclassical eigenvalues 
separately, order by order fl5|j . 

Let us assume that the (n — l) st order approximations 
Wk,n-i to the semiclassical eigenvalues have already been 
obtained and the w k , n are to be calculated. The differ- 
ence between the two subsequent approximations to the 
quantum mechanical response function reads 



Jn(w) = 

k 



m k 



m k 



W - Wk,n + i0 W - Wk,n-1 + i0 
m k Aw k ,n 



(W - W k ,n + iO) 2 



(89) 



with Wk ; 



h(whn+w k ,n-i) and Aw k , n = w k , n -Wk, n -i- 



Integration of (| 

g n (w) = w n 



and multiplication by w n yields 
-m k w n Awk.n 



g n (w)dw = ^ • 



Wh- 



(90) 



which has the functional form of a quantum mechanical 
response function but with residues proportional to the 
order corrections Awk, n to the semiclassical eigenval 



ues. The semiclassical approximation to (90) is obtained 
from the term g n (w) in the periodic orbit sum (^9|) by 
integration and multiplication by it)", i.e., 



Q n (w) = w n J g n (w)dw 



~po 



0[-) . (91) 



We can now Fourier transform both ( |90| ) and (^Tj), and 
obtain (n > 1) 



C n {s) 



2tt 



+ OC 



G n {w)e- iws dw 



iJ2m k (wk) n Awk, n e- iWkS (92) 



= -z£— Aii5(s-s po ). 



(93) 



po 



Equations (|9^) and (^3|) imply that the h expansion of the 
semiclassical eigenvalues can be obtained, order by order, 
by harmonic inversion (h.i.) of the periodic orbit signal in 
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Fig. 7. Correction terms to the semiclassical eigenvalues of 
the circle billiard. Squares: corrections Aw k ,i ~ Wfc.i — Wfe.o 
between first and zeroth order approximations (times multi- 
plicities) obtained by harmonic inversion. Crosses: differences 
ui cx — uiebk between exact quantum and EBK eigenvalues 
(times multiplicities) . 



( P3[ ) to the functional form of (|9^) . [In practice, we will 
again convo lute both expressions with a Gaussian function 
(cf. Section 5.1 ) in order to regularize the 6 functions in 
(|93"|).] The frequencies w k of the periodic orbit signal ( |9^ ) 
are the semiclassical eigenvalues, averaged over different 
orders of h. Note that the accuracy of these values does 
not necessarily increase with increasing order n. We indi- 
cate this in (|92) by omitting the index n at the eigenvalues 
Wk ■ Our numerical calculations for the first order h correc- 
tions show that, in practice, the frequencies Wk we obtain 
are approximately equal to the zeroth order h eigenvalues 
rather than the exact average between zeroth and first 
oder eigenvalues. The corrections Awk, n to the eigenval- 
ues are not obtained from the frequencies, but from the 
amplitudes, mk(wk) n Awk, n , of the periodic orbit signal. 

We will now apply the above technique to the circle 
billiard to obtain the first order corrections to the semi- 



classical eigenvalues obtained in Section 5.2.2 . In Section 
|2.2| , we derived the zeroth order amplitudes of the circle 
billiard (cf. Eq. (Ei6|) 



1 



AW 



2 TOm mT 



3 /2 



(94) 



with sm and ttim the action and multiplicity of the orbit, 
respectiv ely. The fir st order amplitudes are given by (cf. 
Sections O and O): 



^ p ° ~ V r 6sin 3 / 2 



-i(§M r 7T-f) 



7 



(95) 



with 7 = TrM v /M r . 

Using these expressions, we have calculated the first 
order corrections Aw ky i to the lowest eigenvalues of the 
circle billiard, by harmonic inversion of periodic orbit sig- 
nals with s max = 200. Part of the spectrum is presented in 
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Fig. 8. Semiclassical errors \wk,o~ w e x| (diamonds) and |wfc,i — 
w e x| (crosses) of zeroth and first order approximations to the 
eigenvalues obtained by harmonic inversion, in units of the 
average level spacing {Aw) av « 4/to. States are labeled by the 
quantum numbers (n,m). 



Figure [?]. The peak heights (squares) are the corrections 
Awk,i = Wk.i — Wk,o times the multiplicities. For com- 
parison, the differences between the exact and the EBK 
eigenvalues at the positions of the EBK eigenvalues are 
also plotted (see the crosses in Fig. [?]). Both spectra are 
in excellent agreement. The small deviations of the peak 
heights arise from second or higher order corrections to 
the eigenvalues. 

An appropriate measure for the accuracy of semiclas- 
sical eigenvalues is the deviation from the exact quan- 
tum eigenvalues in units of the average level spacings, 
(Aw) av = l/p(w). Figure || presents the semiclassical er- 
ror in units of the average level spacings {Aw) &v rs A/w 
for the zeroth order (diamonds) and first order (crosses) 
approximations to the eigenvalues. In zeroth order approx- 
imation the semiclassical error for the low lying states is 
about 3 to 10 percent of the mean level spacing. This error 
is reduced in the first order approximation by at least one 
order of magnitude for the least semiclassical states with 
radial quantum number n = 0. The accuracy of states with 
n > 1 is improved by two or more orders of magnitude. 



5.4 Reduction of required signal length via harmonic 

inversion of cross-correlated periodic orbit sums 

As described in the sections above, the harmonic inver- 
sion method is able to extract quantum mechanical eigen- 
values from the semiclassical periodic orbit sum includ- 
ing periodic orbits up to a finite action s max - This means 
that in practice, although the periodic orbit sum does not 
converge, the eigenvalues can be obtained from a finite 
set of periodic orbits. The required signal length for har- 
monic inversion depends on the mean density of states, 
i.e., s max w Anp(w) (cf. Depending on the mean 



density of states, the action s max up to which the periodic 
orbits have to be known may therefore be large. Due to 
the rapid proliferation of the number of periodic orbits 
with growing action, the efficiency and practicability of 
the procedure depends sensitively on the signal length re- 
quired. This is especially the case when the periodic orbits 
have to be found numerically. 

The quantization method can be improved with the 
help of cross-correlated periodic orbit sums. The extended 
response functions weighted with prod ucts of diagonal ma- 
trix elements discussed in Section 2.3, in combination with 
the method for harmonic inversion of cross-correlation 



functions presented in Section 3.2, can be used to signif- 



icantly reduce the signal length required for the periodic 

This technique is particularly 
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helpful for chaotic systems, where the periodic orbits must 
be found numerically and where the number of periodic 
orbits grows exponentially with their action. However, for 
regular systems the number of orbits which have to be in- 
cluded can also be significantly reduced as will be demon- 
strated for the circle billiard. 

The basic idea is to construct a set of signals where 
each individual signal contains the same frequencies (i.e., 
semiclassical eigenvalues) and the amplitudes are corre- 
lated by obeying the restriction ([53]) . This can be achieved 
with the help of the gene ralized periodic orbit sum (|3~i| ) 



introduced in Section 2.3 



A set of operators A a , a = 1, . . . , N is chosen. Fol- 
lowing the procedure described in Section 5J , the signals 
C aa > (s) are obtained as Fourier transform of the general- 
ized response functions g^i(w), i.e., 



gZ' H - E AoA a , p A a ,^ w , (96) 

po 

C Qa ' (s) ^ ^ ApoA a pA a ' p6(s Spo) j (^7) 



where for integrable system s the me ans A a , v arc defined by 
(|32|). According to Sections [T^ and 5.1, the corresponding 
quantum mechanical signal is given by 



a 



(• s ) = -iy2m k (k\A a \k) (k\A a ,\k)e- 



qm,oa' 



fc 



(98) 



where t he a mplitudes have the required form ([53|). As in 
Section 5.1 the semiclassical eigenvalues Wk are obtained 
by adjusting the periodic orbit signal (97) (after convolu- 
tion with a Gaussian function) to the functional form of 
the cross-correlated quantum signal (p8|) with the impor- 
tant difference that we now apply the extension of har- 
monic inversion to cross-correlation functions (see Section 
0). 

For the circle billiard, the mean density of states - with 
all multiplicities taken as one - is given by p(w) = w/A. 
According to (fl3|), the signal length required for a single 
signal to resolve the frequencies in a given interval around 
w is therefore approximately given by 



Anp(w) 



2S 



H 



(99) 
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Table 2. Nearly degenerate eigenvalues of the circle billiard, 
obtained by harmonic inversion of a 2 x 2 cross-correlated signal 
of length s max = 150. The denotations are the same as in 
Table hi. The nearly degenerate eigenvalues are now resolved, 
which for a single signal would have required a signal length of 

Smax ~ 500. 



Table 3. Maximum frequencies u> max up to which the spec- 
trum could be resolved using a N x N cross-correlated signal. 
The b a are the operators or functions of operators used to 
build the signal (see text), r: distance from center, L: angular 
momentum in units of hw. 



n 


m 




WEBK 






m hi 


1 


4 


11.064709 


11.048664 


2 


11.048664 


2.0665 





7 


11.086370 


11.049268 


2 


11.049295 


1.9315 


3 


1 


13.323692 


13.314197 


2 


13.314205 


1.9987 





9 


13.354300 


13.315852 


2 


13.315839 


2.0016 



N 


^ max 


1 


45 


2 


65 


3 


90 


4 


120 


5 


130 



1 

1, (r) 

1, <r>, L 2 

1, (r), L 2 , e 

1, <r>, L 2 , l/(r) 2 



((r}-0.7) 2 /3 

(L-l) 2 /10 



10 



where Sh = 27rp is the Heisenberg period (which is action 
instead of time for scaling systems). By using an N x N set 
of signals, it should be possible to extract about the same 
number of semiclassical eigenvalues from a reduced signal 
length s max <C 2Sh , or, vice versa, if the signal length is 
held constant, the resolution and therefore the number of 
converged eigenvalues should significantly increase. 

To demonstrate the power of the cross-correlation tech- 
nique, we first analyze a 2 x 2 cross-correlated periodic 
orbit signal of the circle billiard with A\ = 1 the identity 
operato r and Ai — r. For comparison with the results in 
Section 5.2.2 we choose the same signal length s max = 150. 
By contrast to the eigenvalues in Table p] obtained from 
the one-dimensional signal the nearly degenerate states 
around w ss 11.05 and w rs 13.3 are now resolved as can 
be seen in Table g. Note that a signal length s max ss 500 
is required to resolve these states without application of 
the cross-correlation technique. 

As in all other calculations concerning cross-correlated 
signals, the results were improved by not making a sharp 
cut at the accumulation points but by damping the ampli- 
tudes of the orbits near these points. With the same cut-off 
criterion at the accumulation points, the total number of 
orbits in the signal with s max = 150 was about 10 times 
smaller than in the signal with s max = 500 . This means, 
we could reduce the required number of orbits by one order 
of magnitude. For chaotic systems, where the number of 
orbits grows more rapidly (exponentially) with the maxi- 
mum action, the improvement in the required number of 
orbits may even be better. 

We now investigate the number of eigenvalues which do 
converge for fixed signal length but different sets and di- 
mension of the cross-correlation matrix. Indeed, the high- 
est eigenvalue w max which can be resolved increases sig- 
nificantly when the cross-correlation technique is applied. 
However, the detailed results depend on the operators cho- 
sen. Furthermore, with increasing dimension of the ma- 
trix, the range in which the transition from resolved to 
unresolved eigenvalues takes place becomes broader, and 
the amplitudes in this region become less well converged. 
Rough estimates of w max for various sets of operators and 
fixed signal length s max = 150 are given in Table |. For 
some of the signals, we used the extension of the trace for- 
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Fig. 9. Eigenvalues of the circle billiard in the dense part of the 
spectrum obtained from a 5 x 5 cross-correlated signal of length 
s max = 150 = 0.735iSh with Sh the Heisenberg period. With a 
single signal the required signal length would be s max w 2Sh- 
Squares and sticks: EBK eigenvalues. Crosses: results of har- 
monic inversion. The peak heights give the semiclassical error 
|webk — i0ex| (squares) and the difference between harmonic 
inversion results and EBK eigenvalues |wm — uiebk| (crosses) in 
units of the mean level spacing {Aw) av w 4/ui. The error of the 
harmonic inversion procedure is about an order of magnitude 
smaller than the semiclassical error. 



mula to functions of matrix elements, Eqs. (|35| ) and (|36|) . 
The improvement achieved by increasing the dimension of 
the matrix by one is most distinct for very small N; for 
N > 5, the improvement is only small. This suggests that, 
at given signal length and frequency range, the matrix di- 
mension should not be chosen too large, i.e., there exists 
an optimal matrix dimension, which at constant signal 
length becomes larger with increasing eigenvalues w. 

With a 5 x 5 signal of length s max = 150, eigenvalues 
up to the region w « 130 can be resolved. The results are 
presented in Figure g. There are two points which should 
be emphasized: The first point is that, even in this dense 
part of the spectrum, the error of the method is still by 
about one order of magnitude smaller than the semiclassi- 
cal error, which is illustrated in Fig. ||by the peak heights. 
The squares and crosses mark the semiclassical error 
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I^ebk — w cx \ and the numerical error |u>hi — webk| of the 
harmonic inversion procedure in units of the mean level 
spacing (z4u>) av w A/w. The second point concerns the sig- 
nal length compared to the Heisenberg action Sh = 2np. 
For w = 130, one obtains Sh ~ 204.2. A one-dimensional 
signal would have required a signal length s max s» 2Sh- 
With the cross-correlation technique, we calculated the 
eigenvalues from a signal of length s max = 150 « 0.735S^. 
This is about the same signal length as required by the 
semiclassical quantization method of Berry and Keating 
]l0| , which, however, only works for ergodic systems. 

In summary, our results demonstrate that by ana- 
lyzing cross-correlated signals instead of a single signal, 
the required signal length can indeed be significantly re- 
duced. Clearly, the signal length cannot be made arbitrar- 
ily small, and the method is restricted to small dimensions 
of the cross-correlation matrix. However, the number of 
orbits which have to be included can be very much re- 
duced. Another advantage of the method is that not only 
the frequencies and the multiplicities but also the diago- 
nal matrix elements of the chosen operators are obtained 
by one single calculation. 



5.4.1 Including higher order h corrections 

In the cases discussed so far, we have constructed the 
cross-correlated signal by including different operators and 
making use of Eq. (53). By this procedure, we could ob- 
tain the semiclassical eigenvalues from a signal of reduced 
length or improve the resolution of the spectrum at con- 
stant signal length, while simultaneously obtaining the di- 
agonal matrix elements of the operators. We can now even 
go one step further and include higher H correctio ns i n the 
signal. Here we make use of the results of Sec tion 5.3. The 
first order correction term, which in Section 5.3 was har- 



Table 4. Zeroth (wk,a) and first (Wk,i) order semiclassical ap- 
proximations to the eigenvalues of the circle billiard, obtained 
simultaneously by harmonic inversion of a 3 x 3 cross-correlated 
signal of length s max = 150. The nearly degenerate eigenvalues 
atw« 11.0 and w ~ 13.3 are well resolved. 



monically inverted as a single signal, is now included as 
part of a cross-correlated signal. This procedure combines 
all the techniques developed in the previous sections. 

Formally, the frequencies in the zeroth and first order 
h parts of the cross-correlated signal are not exactly the 
same [see the denominators in Eqs. (75h and (90)], how- 
ever, as already mentioned in Section 5.3, the numerically 
obtained values for the frequencies Wk in ( |92| ) are equal 
to the lowest order h eigenvalues rather than the exact 
average of the zeroth and first order eigenvalues. In prac- 
tice, the cross-correlated signal is therefore in fact of the 
form (|52}). We can now, on the one hand, improve the res- 
olution of the spectrum, and, on the other hand, obtain 
semiclassical matrix elements and the first order correc- 
tions to the eigenvalues with the same high resolution by 
one single harmonic inversion of a cross-correlated signal. 

As an example, we built a 3 x 3 signal for the circle bil- 
liard from the first order correction term given by ( p5| ) and 
the operators A\ = 1 (identity) and A2 = f. Again, we 
chose a signal length of s max = 150. By harmonic inversion 
of the cross-correlated signal, we obtained simultaneously 
the semiclassical eigenvalues, their first order order cor- 
rections, and the semiclassical matrix elements of the op- 
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WEBK 


Wk,0 


w k ,i 


^ cx 








2.356194 


2.356230 


2.409288 


2.404826 





1 


3.794440 


3.794440 


3.834267 


3.831706 





2 


5.100386 


5.100382 


5.138118 


5.135622 


1 





5.497787 


5.497816 


5.520550 


5.520078 





3 


6.345186 


6.345182 


6.382709 


6.380162 


1 


1 


6.997002 


6.997006 


7.015881 


7.015587 





4 


7.553060 


7.553055 


7.590990 


7.588342 


1 


2 


8.400144 


8.400145 


8.417503 


8.417244 


2 





8.639380 


8.639394 


8.653878 


8.653728 





5 


8.735670 


8.735672 


8.774213 


8.771484 


1 


3 


9.744628 


9.744627 


9.761274 


9.761023 





6 


9.899671 


9.899660 


9.938954 


9.936110 


2 


1 


10.160928 


10.160949 


10.173568 


10.173468 


1 


4 


11.048664 


11.048635 


11.063791 


11.064709 





7 


11.049268 


11.049228 


11.087943 


11.086370 


2 


2 


11.608251 


11.608254 


11.619919 


11.619841 


3 





11.780972 


11.780993 


11.791599 


11.791534 





8 


12.187316 


12.187302 


12.228037 


12.225092 


1 


5 


12.322723 


12.322721 


12.338847 


12.338604 


2 


3 


13.004166 


13.004167 


13.015272 


13.015201 


3 


1 


13.314197 


13.314192 


13.323418 


13.323692 





9 


13.315852 


13.315782 


13.356645 


13.354300 


1 


6 


13.573465 


13.573464 


13.589544 


13.589290 


2 


4 


14.361846 


14.361846 


14.372606 


14.372537 





10 


14.436391 


14.436375 


14.478531 


14.475501 


3 


2 


14.787105 


14.787091 


14.795970 


14.795952 


1 


7 


14.805435 


14.805457 


14.821595 


14.821269 


4 





14.922565 


14.922572 


14.930938 


14.930918 





11 


15.550089 


15.550084 


15.593060 


15.589848 


2 


5 


15.689703 


15.689701 


15.700239 


15.700174 


1 


8 


16.021889 


16.021888 


16.038034 


16.037774 


3 


3 


16.215041 


16.215047 


16.223499 


16.223466 


4 


1 


16.462981 


16.462982 


16.470648 


16.470630 





12 


16.657857 


16.657846 


16.701442 


16.698250 


2 


6 


16.993489 


16.993486 


17.003884 


17.003820 


1 


9 


17.225257 


17.225252 


17.241482 


17.241220 


3 


4 


17.607830 


17.607831 


17.615994 


17.615966 





13 


17.760424 


17.760386 


17.804708 


17.801435 


4 


2 


17.952638 


17.952662 


17.959859 


17.959819 


5 





18.064158 


18.064201 


18.071125 


18.071064 



erator f. The results for the zeroth order approximations 
Wk,o to the eigenvalues and the first order approximations 
Wk,\ = Wk,Q-\-Awk,i are presented in Table For compar- 
ison the exact and the EBK eigenvalues are also given. As 
for the results presented in Table ||, we were able to resolve 
the nearly degenerate states in the zeroth order approx- 
imation, which for a single signal would have required a 
signal length of s max 500. Moreover, in contrast to the 
results of Section |5.3| , we could now also resolve the first 
order approximations to the nearly degenerate states. 
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6 Conclusion 

The harmonic inversion method has been introduced as a 
powerful tool for the calculation of quantum eigenvalues 
from periodic orbit sums as well as for the high resolution 
analysis of quantum spectra in terms of classical periodic 
orbits. We have demonstrated that this method, which has 
already successfully been applied to classically chaotic sys- 
tems, yields excellent results for regular systems as well. 
Harmonic inversion has thus been shown to be a universal 
method, which, in contrast to other high resolution meth- 
ods, does not depend on special properties of the system 
such as ergodicity or the existence of a symbolic code. 

With the harmonic inversion method, we are able to 
calculate the contributions of the classical periodic orbits 
to the trace formula from the quantum eigenvalues with 
high precision and high resolution. By analyzing the dif- 
ference spectrum between exact and semiclassical eigen- 
values, we could determine higher order ft corrections to 
the periodic orbit sum of the circle billiard. 

Up to now, no theory for the ft corrections to the Berry- 
Tabor formula for regular systems has been developed. We 
have numerically found an expression for the first order ft 
corrections to the Berry- Tabor formula by harmonic inver- 
sion of the difference spectrum. The same expression can 
be derived analytically by using Vattay's and Rosenqvist's 
method for chaotic systems and introducing some reason- 
able ad-hoc assumption for the circle billiard. As this is 
clearly not a strict derivation, it is an interesting task for 
the future to develop a general theory for the higher order 
ft corrections to the trace formula for regular systems. 

In addition to calculating semiclassical eigenvalues 
from the usual periodic orbit sum, we have demonstrated 
how further information can be extracted from the pa- 
rameters of the classical orbits by applying the harmonic 
inversion technique to different extensions of the trace for- 
mula. Using a generalized trace formula including an arbi- 
trary operator, we have shown that the method also allows 
the calculation of semiclassical diagonal matrix elements 
from the parameters of the periodic orbits. Furthermore 
we have demonstrated how higher order ft corrections to 
the semiclassical eigenvalues can be obtained by harmonic 
inversion of correction terms to the periodic orbit sums. 
For the case of the circle billiard, we found that, including 
the first order correction, the accuracy of the semiclassical 
eigenvalues compared to the exact quantum eigenvalues 
could be improved by one or more orders of magnitude. 

Although by harmonic inversion the quantum eigenval- 
ues can be calculated from a semiclassical signal of finite 
length, i.e., from a finite set of periodic orbits, the number 
of orbits which have to be included may still be large. We 
have demonstrated that by a generalization of the har- 
monic inversion method to cross-correlation functions the 
required signal length may be significantly reduced, even 
below the Heisenberg time. Because of the rapid prolifer- 
ation of periodic orbits with growing period, this means 
that the number of orbits which have to be included may 
be reduced by about one to several orders of magnitude. 



A Calculation of the first order h correction 
terms to the semiclassical trace formula 

The calculation of higher order ft correcti ons to the semi- 
classical eigenvalues introduced in Section 2.4 requires the 
knowledge of the nth order amplitudes A.p$ in the peri- 
odic orbit sum (j39|). In this Appendix, we briefly outline 

the derivation of the first order amplitudes A.p~o and the 
application to the circle billiard given in Eq. (|Io|). 

Two different methods for the calculation of higher 
order ft correction terms in chaotic systems have been 
derived by Gaspard and Alonso 31,3^] and Vattay and 
Rosenqvist P3,pj,p5|. The latter method has been spe- 



cialized to two-dimensional chaotic billiards in 35 . Here, 
we follow the approach of Vattay and Rosenqvist. How- 
ever, it is important to note that both methods cannot 
straightforwardly be applied to integrable systems and ad- 
ditional assumptios will be necessary to derive Eq. (fib] ) for 
the circle billiard. A general theory for the calculation of 
higher order ft corrections to the Berry- Tabor formula (^) 
for integrable systems is, to the best of our knowledge, 
still lacking. 

Vattay and Rosenqvist give a quantum generalization 
of Gutzwiller's trace formula based on the path integral 
representation of the quantum propagator. The basic idea 
of their method is to express the global eigenvalue spec- 
trum in terms of local eigenvalues computed in the neigh- 
bourhood of periodic orbits. The energy domain Green 
function G{q 1 q' 1 E) is connected to the spectral determi- 
nant A(E) = II n (E — E n ) , with E n the energy eigenvalues 
or resonances, by 



Tr G(E) = / dqG(q,q,E) 



dE 



In A(E). 



(100) 



The trace of the Green function can be expressed in terms 
of contributions from periodic orbits 



Tr G(E) = Tr G P (E), 



(101) 



p.o. 



with the local traces connected to the local spectral de- 
terminants by 



Tr G P (E) = — In A p (E). 



(102) 



The trace of the Green function can therefore be calcu- 
lated by solving the local Schrodinger equation around 
each periodic orbit, which yields the local eigenspectra. 
To obtain the local eigenspectra, the ansatz 



(103) 



is inserted into the Schrodinger equation, yielding the fol- 
lowing differential equations for <P and S. 

d t S +^(\7S) 2 + U = (104) 

<9 t + V<£VS'+i<Z>AS-yZ\<Z> = O, (105) 
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where U is the potential entering the Schrodinger equa- 
tion. 

The spectral determinant can be calculated from the 
local eigenvalues of the amplitudes <!>. For arbitrary energy 
E, the amplitudes <P l p of the local eigenfunctions fulfill the 
equation 

-P p (t + T p ) = R l p (E)<P p (t), (106) 

wher e T v is the period of the classical orbit. Using 
Eq. (102), the trace formula can be expressed in terms 
of the eigenvalues R l p (E): 

1 ( d\nRL(E)\ 

oo 



(107) 



r=l 



This is the quantum generalization of Gutzwiller's trace 
formula and holds exactly. 

The amplitudes and their eigenvalues are now ex- 
panded in powers of h: 



(108) 

(109) 
(110) 




exp(Cn(l+2-q 



The terms yield the Gutzwiller trace formula as ze- 
roth order approximation, while the terms with m > 
give H corrections. 



To solve Eqs. (104) and ( Jl 05| ) in different order of h, 
the Schrodinger equation and the functions ^( m ) and 5* 
are Taylor expanded around the periodic orbit, 

S(q.i) = E^[ s «(*)( c l -%(*))" (HI) 

^)(q,t) = ^^L M (*)(q-qp(*)) n , (H2) 

resulting in a set of differential equations for the different 
orders of the Taylor expansions and different orders in h. 
In one dimension, these equations read explicitly: 



1 - 

*n+l?+^E 



jrys n _ fc+ is fc+ i+u„ = 0, (113) 



2^ (n-k)lk 

k=0 y ' 

where u n are the coefficients of the Taylor expanded po- 
tential, and 



(m+l) 



dm+1) • 
°n+l 1 



E 



fe=0 



k)\k\ 



,(m+l) 



1 







Am) 
°n+2 



(114) 



This set of differential equations can be solved iteratively. 
The Z-th eigenfunction is characterized by the condition 



dm) 



for n < I. The different orders of h are connected 



by the last term in (114). Starting from zeroth order h 
and the lowest nonvanishing order of the Taylor expan- 
sion, the functions can be determined order by order. For 
higher dimensional systems, the coefficient matrices obey 
similar equations, and the structure of the set of equations 
remains the same. 

To obtain the first order h correction to the Gutzwiller 
trace formula, one has to calculate the quantities Cj . To 
obtain the se quantities one has to solve the set of equa- 
tions (114) up to the lowest nonvanishing first order h 

coefficient function <j)\ , respectively. As (jhi" 1 ^ = for 
n < I, this involves solving the equations for S2, S3 and S4, 
and for the zeroth order h coefficient functions 0; < "°\ (frl+l 
and 0;+2- II tne initial conditions are set to be <P L ^°\o) = 1 



and ^ m \o) = for m > 0, the correction term C l 
then given by the relation 



(i) 



exp(Cf 0) ) 



(115) 



whic h follows from the h expansion of the eigenequation 

An explicit recipe for the calculation of the first h cor- 
rection for two-dimensional chaotic billiards is given in 
[E35J. For billiards, the potential U in the Schrodinger equa- 
tion equals zero between two bounces at the hard wall. 
The functions S and <P now have to be Taylor expanded 
in two dimensions: 

S(x, y, t) = S a + S x Ax + S y Ay 

+ 1 (S X 2 (Ax) 2 + 2S xy AxAy + S y *{Ay) 2 ) 
+ ... (116) 

and similarly for If the coordinate system is chosen in 
such a way that x is along the periodic orbit and y is 
perpendicular to the orbit, derivatives with respect to x 
can be expressed in terms of the derivatives with respect 
to y using the stationarity conditions 



S 



J x Tly7> 

Sx 



(117) 



The quantity S x is equal to the classical momentum of the 
particle. 

For the free motion between the collisions with the 



wall, the s et of differential equations correspoding to ( |113| ) 
and (114) then reduces to a set of equations involving 



only derivatives with respect to y. These equations can be 
solved analytically, with the general solution still contain- 
ing free parameters. Setting S x = 1, the first coefficient 
functions of the Taylor expanded phase are given by: 



Syy(t) 



1 



t + t Q 



(118) 
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Syyy(t) 



,(*) 



A 



3 



B 



3A 2 



(t + t ) 3 (t + t y {t + t y 



(119) 
(120) 



where to, A and B are free parameters. For given I, the 
first nonvanishing coefficients of the amplitude read 



AO) 



(t)=E 



to 



t + ta 
E 



1+1/2 



(t + t y+ 3 / 2 

E 

(t+t y+ 5 / 2 



C + (l + l) [ 



D 



A t 



1+1/2 



2 (t + t ) 



(121) 
(122) 



t + t Q 



AH l +1/2 ^ 



2(t + t ) 2 

i(Z + 2) 2 (Z + l) 2 + |(Z + 2)(Z + 1)(| + 

(123) 

Again, D and E are free parameters. 

At the collisions with the hard wall, the phase and 
amplitude have to obey the bouncing conditions 

S(x,y,t_ ) = S(x,y,t +0 ) + iir (124) 
$(x,y,t- ) = $(x,y,t +0 ), (125) 

from which the bouncing conditions for the coefficients of 
the Taylor expanded functions S and <P can be derived. 
While the general solutions between the bounces are valid 
for all billiards, the bouncing conditions in their Taylor 
expanded form depend explicitely on the shape of the hard 
wall. 

An additional condition which the solutions S and <P 
have to obey is periodicity along the orbit. With every 
traversal the phase gains the same constant contributions 
at the collisions with the wall. The derivatives of the phase 
are periodic. The amplitude collects the same factor with 
each traversal, which means that all Taylor coefficients of 
the amplitude are periodic apart from a constant factor. 
These conditions together with the bouncing conditions 
determine the values of the free constants in the general 
solutions between the collisions. 

The solutions can in general be found numerically, by 
choosing suitable initial conditions and following the evo- 
lution of the phase and amplitude functions along the or- 
bit. After several iterations around the orbit, the param- 
eters should converge against their periodic solution. The 
correction terms C, are then given by the integral 



(i) _ 



/ 



yl x Z 



HO) 



ctt, 



(126) 



which can be computed explicitly from the solutions found 
above. 

As already explained, this method is designed for 
chaotic systems, as its derivation is based on the assump- 
tion that the periodic orbits are isolated. Nevertheless, we 
obtain reasonable results when applying the method to the 
circle billiard, taking one periodic orbit from each rational 
torus and introducing some additional assumptions. 

Because of the symmetry of the orbits, we can assume 
that every side of the orbit contributes in the same manner 
to the h correction term for the whole orbit. This means, if 
we reset t = at the start of e ach s ide, the free parameters 
in the general solutions (118) to (123) must be the same 
for each side, apart from the parameter E, which collects 
the same factor during every collison with the wall. With 
these assumptions, the differential equations can be solved 
analytically. 

However, i t tur ns out tha t the bouncing conditions re- 
sulting from (124) and (125) are not sufficient to deter- 
mine all free parameters, as some of the conditions are 
automatically fulfilled. We need additional conditions for 
the parameters. These can be obtained from the rotational 
symmetry of the system: Because of this symmetry, we 
can assume that the amplitude of the wave function does 
not depend on the polar angle ip. The same holds for all 
derivatives of the amplitude with respect to the radius r. 
For the zeroth order h amplitudes, expressed in polar co- 
ordinates (r, ip) , this gives us the additional conditions we 
need: 

If we further assume that the phase separates in polar 
coordinates 

S(r,<p)=S r (r) + S v (<p), (128) 

which implies that all mixed derivatives vanish, it turns 
out that we do not need the bouncing conditions at all. 
All parameters can be determined from the symmetry of 
the system, and the bouncing conditions are automatically 
fulfilled. We c onsidered only the case I — 0, for which we 
used Eq. (128) together with the conditions 



— <P m = 0, 

dtp 



dip dr 



= 0. 



(129) 



Ou r fina l results for the constant parameters in Eqs. ( Jl 18[ ) 
to (p3l) are 

(130) 
(131) 
(132) 
(133) 



<o = — sin 7 
A = — cos 7 
B = 
C = 

D = --sin 1/2 7 
2 ' 



(134) 



with 7 as defined in Section 2.4 (see Fig. [j]). The radius 
of the billiar d w as taken to be R = 1. Inserting these 
solutions in (|126|) finally leads to 



(135) 



= M r 



3 sin 7 



6 sin 3 7 
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where M r is the number of sides of the orbit. 

The first order ampli tude s are obtained b y in- 

serting the h expansion flllOj ) in the trace formul a (Jl07|) 
and comparing the result with the h expansion (|39[). In 
the units we have used here (radius R = 1 and momen- 
tum hk — 1), the scaling parameter w is equal to H. If 
we use only the I = contributions and assume that the 
terms exp(C , Q°' ) ) are equal to the amplitudes given by the 
Berry- Tabor formula, we finally end up with the expres- 
sion (Eoj). Although we cannot strictly justify the l ast s tep, 



pro- 



our analysis of the quantum spectrum in Section 4.2 
vides strong numerical evidence that Eq. ( p| ) is correct. 
It will be an interesting task for the future to develop a 
general theory for the h correction terms of integrable sys- 
tems and thus to provide a more rigorous mathematical 
proof of Eq. (ph. 
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